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I. INTRODUCTION 

Ising model in a transverse field has been widely ex- 
amined in statistical mechanics and condensed matter 
physics since the pioneering work of de Gennes [l] who in- 
troduced it as a pseudo spin model for hydrogen-bonded 
ferroelectrics such as the KH2PO4 type. Following stud- 
ies has been predicated that this semi-quantum mechan- 
ical model can be successfully applied to a variety of 
physical systems such as DyVOi, TbVO^ Q and some 
real magnetic materials [3|. From the theoretical point 
of view, transverse Ising model (TIM) has been investi- 
gated by a variety of techniques such as renormalization 
group method (RG) [|, effective field theory (EFT) 0- 
.7], cluster variation method (CVM) Q, mean field the- 
ory (MFT) pair approximation (PA) [13] and Monte 
Carlo simulations (MC) [11]. In the previous works men- 
tioned above, the authors focused their attention on the 
behavior of tricritical points but have not considered the 
effect of the crystal field (i.e. single ion anisotropy) in 
Hamiltonian describing the system. 

However, there are a few studies in the literature that 
include the crystal field as well as the transverse field 
interactions. Recently, the effect of both the transverse 
field and the crystal field on the spin-S" Ising model with 
spins of magnitude S=l have been studied and it is shown 
that TIM model presents a rich variety of critical phe- 
nomena. For example, Jiang et al. [12| has studied the 
spin-1 TIM on a honeycomb lattice with a longitudinal 
crystal field and discussed the existence of a tricritical 
point at which the phase transition change from second- 
order to first order. By using the EFT with a probability 
distribution technique, Htoutou et al. [l3| have inves- 
tigated the influence of the crystal field on the phase 
diagrams of a site diluted spin-1 TIM on a square lat- 
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tice. Similarly, Jiang have studied a bond diluted spin-1 
TIM with crystal field interaction for a honeycomb lat- 
tice within the framework of the EFT with correlations 
fl4j . In these studies, the authors have reported the ob- 
servation of a reentrant behavior on the system. Further- 
more, in a series of papers Htoutou et al. [HI, [HI have 
discussed the dependence of the behavior of the order 
parameters on the transverse and crystal fields, but they 
have restricted themselves on the second order transition 
properties. The effect of a longitudinal crystal field on 
the phase transitions in spin-3/2 and spin-2 transverse 
Ising model has been also examined for both honeycomb 
and square lattices by using the EFT with correlations 
[l7l fig . More recently, within the basis of EFT and 
MFT, Miao et al. [l9| have studied the phase diagrams 
of a spin-1 transverse Ising model for a honeycomb lat- 
tice. They have obtained the first-order transition lines 
by comparing the Gibbs free energy. 

An ordinary EFT approximation includes spin-spin 
correlations resulting from the usage of the Van der Waer- 
den identities and provides results that are much superior 
to those obtained within the traditional MFT. However, 
the EFT approximations mentioned above are not suffi- 
cient enough to improve the results much. The reason 
may be due to the usage of a decoupling approxima- 
tion that neglects the correlations between different spins 
that emerge when expanding the identities. According 
to us, the first-order transition lines obtained in Refs. 
[l2T - [l6j and [l9| are incomplete, because the spin corre- 
lation functions such as (SqSi), (SqSiS%) etc. can not 
be determined by using any decoupling approximation 
or in other words, by neglecting the correlations between 
different spins. In order to overcome this point, we pro- 
posed the IEFT approximation that takes into account 
the correlations between different spins in the cluster of 
considered lattice (q — 3) [20-22]. Namely, the hallmark 
of the IEFT is to consider the correlations between differ- 
ent spins that emerge when expanding the identities and 
this method is superior to conventional mean field the- 



2 



ory and the other EFT approximations in the literature. 
Therefore, it is expected that the calculation results will 
be more accurate. 

As far as we know, there is not such a study which in- 
cludes the longitudinal component of the magnetic field 
in addition to the crystal field and transverse magnetic 
field on the Hamiltonian. Thus, in the present work, we 
intended to investigate the thermal and magnetic proper- 
ties of spin-1 TIM with crystal field under a longitudinal 
magnetic field on a honeycomb lattice within the frame- 
work of the IEFT. For this purpose, we investigated the 
proper phase diagrams, especially the first-order transi- 
tion lines that include reentrant phase transition regions 
and we improved the results in Refs.[l2|,[l5j]. We gave the 
numerical results for the behavior of the order parame- 
ters when the system undergoes a first or second order 
transition at a finite temperature. In addition, it would 
be interesting to see how the thermodynamic properties 
like entropy S which has not been calculated before and 
Helmholtz free energy F are effected by the crystal field 
as well as transverse and longitudinal magnetic fields. 
Hence, the numerical results are presented and compared 
with the literature. 

The layout of this paper is as follows. In section [Til 
we briefly present the formulations of the IEFT. The re- 
sults and discussions are presented in section ITTTl Finally, 
section HVl contains our conclusions. 



II. FORMULATION 

As our model we consider a two dimensional lattice 
which has N identical spins arranged. We define a clus- 
ter on the lattice which consists a central spin labeled 
So, and q perimeter spins being the nearest-neighbors of 
the central spin. The cluster consists of (q + 1) spins 
being independent from the value of S. The nearest- 
neighbor spins are in an effective field produced by the 
outer spins, which can be determined by the condition 
that the thermal average of the central spin is equal to 
that of its nearest-neighbor spins. The Hamiltonian of 
the spin-1 transverse model with crystal field in a longi- 
tudinal magnetic field is given by 

<i,j> i i i 

where Sf and Sf denote the z and x components of the 
spin operator, respectively. The first summation in equa- 
tion ([1]) is over the nearest-neighbor pairs of spins and 
the operator Sf takes the values Sf = 0, ±1. J, D, Q 
and h terms stand for the exchange interaction, single- 
ion anisotropy (i.e. crystal field) and transverse and lon- 
gitudinal magnetic fields, respectively. 

At first, we start constructing the mathematical back- 
ground of our model by using the approximated spin 
correlation identities introduced by SaBarreto, Fittipaldi 



and Zeks |2J 

({fi}S?) = 



{fi} 



TriS?exp(-PHi) 
Tn exp(-(3Hi) 



({fi}(S?) 2 ) = ({f i } 



TniSffeigyj-pHi) 
Tr i exp(-/3i7 l ) 
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where ft = 1/ksT and a = z or x. 

In order to apply the differential operator technique, 
we should separate the Hamiltonian ([1} into two parts as 
H = Hi + H . Here, one part denoted by Hi includes 
all contributions associated with the site i, and the other 
part H does not depend on the site i. At this point, one 
should notice that Hi and H do not commute with each 
other. We can write —Hi as 



H = E,sf + d (sty + nsf + hSf 



(4) 



where Ei — JJ2j &j ^ s the local field on the site i. If we 
use the matrix representations of the operators Sf and 
Sf for the spin-1 system then we can obtain the matrix 
form of equation 



Ei + D + h Q/V2 

//, - ! n/V2 o n/Vz 

Q/V2 -Ei + D-h 

In order to proceed further, we have to diagonalize 
matrix in equation ([5]). The three eigenvalues are 
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D + 18Eih + 9fr 



D 2 + 6Eih + 3h 2 , 



and the eigenvectors (p k of —Hi corresponding to the 
eigenvalues in equation (JSJ are calculated as follows 
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Hereafter, we apply the differential operator technique in 
equations © and © with {/j} = 1. From equation © 
we obtain the following spin correlations for the thermal 
average of a central spin for honeycomb lattice (q = 3) 
as 

{S$) = / JJ [1 + 5Jsinh(JV) + (5 , |) 2 {cosh(JV) - 1}] \ 
xF(x)\ x=0 , (8) 

(5£) = / JJ [1 + 5?sinh(JV) + (5|) 2 {cosh(JV) - 1}] \ 

xH(x)\ x=0 . (9) 

By expanding the right-hand sides of equations © and 
© we get the longitudinal and transverse spin correla- 
tions as 

m z = (5 z ) = l + 3fci(Si) + 3(Ii - lo)(S!) + 3/ 2 <5i5 2 ) 

+6(fc 2 - fciX^iSf > + 3(; - 2h + h)(s!si) 

+k 3 {S 1 S 2 S 3 ) + 3(Z 4 - ; 2 )(^5 2 5 3 2 ) 

+3(fei -2k 2 + k i ){S l S 2 2 Sl) 

+ (-/ + 3/i - 3/ 3 + ^XS 2 ^ 2 ), (10) 

- (5ff) = po + 3 Cl (^) + 3(pi -ft>)(5 2 ) + 3p 2 <5 x 5 2 ) 
+6(c 2 - cO^Sf) + 3( Po - 2pi +P 3 )(5 2 5 2 2 ) 
+c 3 (5i5 2 5 3 ) + 3(p 4 - P2){S 1 S 2 Sl) 
+3{c 1 -2c 2 + c i )(S 1 S 2 2 S 2 3 ) 
+(-po + 3pi - 3p 3 + P5)(5 2 5 2 5 3 2 ). (11) 

Next, the average value of a perimeter spin in the system 
can be written as follows and it is found as 

mi = (5f ) = (1 + ^sinh(JV) + (5 z ) 2 {cosh( JV) - 1}> 
xF(z + 7 )U = o, (12) 

(S 1 ) = ai (1 - ((^) 2 » + a 2 <5 *) + a 3 ((^) 2 ), (13) 

where 7 = (3 — 1)A is the effective field produced by 
the (q — 1) spins outside the system and A is an un- 
known parameter to be determined self-consistently. In 
the effective-field approximation, the number of indepen- 
dent spin variables describes the considered system. This 
number is given by the relation v — ((S^) 2S ). As an ex- 
ample for the spin-1 system, 25 = 2 which means that 
we have to introduce the additional parameters ((5q) 2 ), 
((5q) 2 ) and ((5|) 2 ) resulting from the usage of the Van 
der Waerden identity for the spin-1 Ising system. With 
the help of equation ([3]) 

q, = ((5 z ) 2 ) 

= / JJ [1 + 5?sinh(JV) + (5|) 2 {cosh( JV) - 1}] \ 
xG(x)U =0 , (14) 



& = ((SoT) 

= ^ JJ [1 + 5|sinh(JV) + (5|) 2 {cosh(JV) - 1}] \ 

y.K{x)\ x=Q . (15) 

Hence, we get the quadrupolar moments by expanding 
the right-hand sides of equations (TH| and (|15[) 



((5 2 ) 2 ) = r + 3n 1 (5 1 )+3(r 1 -r )(5 1 2 )+3r 2 (5 1 5 2 ) 

+6(n 2 - nx^SxSfj + 3(r - 2r x + r 3 )(5 2 5 2 ) 
+n 3 (5i5 2 5 3 ) + 3(r 4 - r 2 )(SiS 2 S 2 } 
+3(ni-2n 2 + n 4 )(5i5^5|) 
+ (-r + 3ri - 3r 3 + r 5 )(S 2 S 2 S 2 ), (16) 



((5^) 2 ) = v + 3fi 1 (S 1 ) + 3(v 1 -v )(S 2 1 ) + 3v 2 (S 1 S 2 ) 

+6(p 2 - Mi)<5i5|) + 3(tj - 2«i + « 3 )(5 2 5 2 ) 
+/j 3 (5i5 2 5 3 ) + 3(u 4 - « 2 )(5i5 2 5 2 ) 
+3( M i -2 Ai2 + / J4 )(5 1 5 2 5 3 2 ) 
+ (-t; + 3«i - 3v 3 + v 5 )(S 2 S 2 S 2 ). (17) 

Corresponding to equation (fT2|) 

((5|) 2 ) = <l+5 z sinh(JV)+(5 2 ) 2 {cosh(JV)-l})G( a; +7), 



(18) 

<5 2 ) = h (1 - ((5 *) 2 )) + b 2 {S* ) + 6 3 ((5 2 ) 2 ). (19) 

Details of calculation through ([BT fT^l) can be found in Ap- 
pendix section. With the help of equations © and © 
the functions F(x), G(x), H(x) and K(x) in equations 
© j @ , (O an d (f!5l) can be calculated numerically from 
the relations 



F(x) 



H(x) 



1 



s=3 



a — y2{<Pn\Sf\<p n ) exp(/3A„), 
En=l eX P (/ 3A «) n=l 



(20) 



s=3 



^(^„|5f|(^ n )exp (fiXn), 



(21) 



s=3 



G(a;) = —5=3 

En=l ex P (0*n) n=l 



X>„|(Sf) 2 |¥>n)exp(/3A n ), 

(22) 

K{x) = — ^ — -^(^|(5r) 2 |^)exp(/3A„). 

E„=iexp (/3A„) n=1 

The internal energy U per site of the system can be ob- 
tained easily from the thermal average of the Hamiltonian 
in equation ©. Thus, the internal energy is given by 



NJ 



(5o50 + D((S^) + Q(SS) + h{S%), (24) 
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where the correlation functions (So Si), ((Sq) 2 ), (Sq) and 
(5q) are obtained from equation (|A7|) . With the use of 
equation ([24]) . the specific heat of the system can be nu- 
merically determined from the relation 

The Helmholtz free energy of a system is defined as 

F = U -TS, (26) 

which, according to the third law, can be written in the 
form |25f 

F = U-tJ^ ^dT', (27) 

where the integral in the second term is entropy of the 
system according to the second law. First, we calculate 
the internal energy per site from equation (|24p . Then 
with the help of equation (|25l) we can carry out numerical 
integration and calculate the entropy and free energy of 
the system. Equations (I0J), (JTT]), ([15]), (fT6]) . (ITT)) and 
(H2J) are fundamental correlation functions of the system. 
When the right-hand sides of equations © , © , (IT41 and 
(j!5[) are expanded, the multispin correlation functions 
can be easily obtained. The simplest approximation, and 
one of the most frequently adopted is to decouple these 
equations according to 

<^(^) 2 ...5f)-(5f)<(5|) 2 )...(^), (28) 

for i j ^ ... ^ I [24j. The main difference of the method 
used in this study from the other approximations in the 
literature emerges in comparison with any decoupling ap- 
proximation (DA) when expanding the right-hand sides 
of equations ©, ©, (O and (JT3J!. In other words, 
one advantage of the approximation method proposed by 
this study is that no uncontrolled decoupling procedure 
is used for the higher-order correlation functions. 

For spin-1 Ising system with q = 3, taking equations 
([TO]). (JID , (fl"3"j) . OH), dHJ) and dinj) as a basis, we derive 
a set of linear equations of the spin correlation functions 
which interact in the system. At this point, we assume 
that (i) the correlations depend only on the distance be- 
tween the spins, (ii) the average values of a central spin 
and its nearest-neighbor spin (it is labeled as the perime- 
ter spin) are equal to each other, and (iii) in the matrix 
representations of spin operator S, the spin-1 system has 
the properties (Sf) 3 = Sf and (Sf) 4 = (Sf) 2 . Thus, the 
number of the set of linear equations obtained for the 
spin-1 Ising system with q = 3 reduces to twenty three 
and the complete set is given in Appendix. 

If equation (|A7|) is written in the form of a 23 x 23 
matrix and solved in terms of the variables Xi[{i — 
1,2, ...,23) {e.g., x\ = (Sq z ),x 2 = (S 1 S ), x 23 = 
((Sq) 2 ))} of the linear equations, all of the spin corre- 
lation functions can be easily determined as functions of 



the temperature, effective field, crystal field and longitu- 
dinal magnetic field as well as transverse magnetic field 
which the other studies in the literature do not include. 
Since the thermal average of the central spin is equal 
to that of its nearest-neighbor spins within the present 
method then the unknown parameter A can be numeri- 
cally determined by the relation 

(S5) = (Si) or Xl = Xi . (29) 

By solving equation (|29[) numerically at a given fixed 
set of Hamiltonian parameters we obtain the parameter 
A. Then we use the numerical values of A to obtain 
the spin correlation functions (Sq), (Sq) (longitudinal 
and transverse magnetizations), ((Sq) 2 ), ((Sq) 2 ) (longi- 
tudinal and transverse quadrupolar moments) and so on, 
which can be found from equation (|A7[) . Note that A = 
is always the root of equation (f2"9"| corresponding to the 
disordered state of the system. The nonzero root of A in 
equation (j2"9"|) corresponds to the long-range ordered state 
of the system. Once the spin correlation functions have 
been evaluated then we can give the numerical results for 
the thermal and magnetic properties of the system. 

III. RESULTS AND DISCUSSIONS 

In this section, we can examine the ferromagnetic prop- 
erties of the spin-1 TIM with crystal field under an ap- 
plied longitudinal magnetic field on a honeycomb lattice 
using the IEFT. For this purpose, we focus our attention 
on the phase diagrams of the system in (ksT c / J — Q/ J) 
and (kBT c /J — D/J) planes and investigate the whole 
phase diagrams by examining the numerical results for 
the thermal and magnetic properties. In order to plot 
the phase diagrams, we assume (Sq) = (Si) and the ef- 
fective field 7 is very small in the vicinity of fcsT c / J and 
solve the set of linear equations in equation (|A7|) nu- 
merically using the self-consistent relation correspond- 
ing to equation (f2"9"|) . In Figs, la and lb, we plot 
the variation of the critical temperature with transverse 
field fi/J and crystal field D/J, respectively. Fig. la 
shows the phase diagram in the (UbTc/ J — CI/ J) plane 
with h/J — and for selected values of D/J, namely 
-1.0,-1.3,-1.381,-1.382 and -1.4. In this figure, we 
can call attention to the signs of an interesting behavior 
known as reentrant phenomena. In other words, when 
the crystal field strength is positive valued, the type of 
the transition in the system is invariably second order 
which is independent from transverse field value. On 
the other hand, if the crystal field value is sufficiently 
negative then we can expect to see two successive phase 
transitions. Solid and dashed lines in Fig. la correspond 
to the second and first order phase transition lines, re- 
spectively. Tricritical end points at which first and sec- 
ond order transition points meet are shown as white cir- 
cles. In our calculations, we realized that one can ob- 
serve reentrant behavior in the system for the values of 
Cl/J < 0.861 and -1.4533 < D/J < -1.0201. For the 
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FIG. 1: Phase diagrams of the spin-1 system with h/J = in 
(a) (k B T c /J-tt/J) y (b) (k B T c /J-D/J) planes. The numbers 
on the curves denote the values of the crystal field D/J and 
transverse field fl/ J , respectively, (c) Transverse field depen- 
dencies of the tricritical temperature ksTt/J and tricritical 
crystal field —D t /J. 



values of D/J < —1.382 the transition lines exhibit a 
bulge which gets smaller as the value oi D/J approaches 
the value of —1.4533 which means that ferromagnetic 
phase region gets narrower. We have also examined the 
phase diagram of the present system in (kgT c / J — D / J) 
plane with h/ J — and for selected values of fl / J such as 
0,0.25,0.5,0.75,0.86, 1.0 and 1.1. The numerical results 



are shown in Fig. lb. Solid and dashed lines in Fig. lb 
correspond to the second and first order phase transition 
lines, respectively. White circles denote tricritical points. 
As we can see from this figure, as the value of transverse 
field f2/J increases starting from zero then the value of 
tricritical point decreases gradually and disappears for 
fl/J> 0.86. If the transverse field value is greater than 
this value then we have only second order transitions in 
the system. These results show that the reentrant phe- 
nomenon originates from the competition between the 
crystal field D/J and transverse field fl/J. Furthermore, 
the variation of the coordinates of the tricritical points 
ksTt/J and D t /J as a function of transverse field CI/ J 
is illustrated in Fig. lc. This figure shows that the tri- 
critical points exist for 1.4333 < —Dt/J < 1.2776 and 
51/ J < 0.861. In addition, value of ksT t /J and the ab- 
solute value of D t / J decreases as the value of transverse 
field increases and the tricritical temperature disappears 
at the critical value of the transverse field fit/ J — 0.861. 
All of the results mentioned above are in a good agree- 
ment with other works M, fl3-fl6ll , but not with Qjl . 



In order to clarify the first-order phase transitions in 
the system we examine the variation of the order pa- 
rameters with temperature such as the longitudinal and 
transverse magnetizations m z and m x together with the 
longitudinal and transverse quadrupolar moments q z and 
q x . The effects of the crystal field on the behavior of the 
order parameters with a selected transverse field value 
SI/ J = 0.5 are shown in Fig. 2. Here, dashed lines cor- 
respond to the solutions of the first order transition. As 
we can see on the left panel in Fig. 2a, for the selected 
values of the crystal field D/J — 1.0,0,-0.5,-1.0 and 
— 1.25 with h/J = as the temperature increases, the 
longitudinal magnetization m z falls rapidly from its sat- 
uration magnetization value at fc^T/J = and decreases 
continuously in the vicinity of the transition tempera- 
ture and vanishes at a critical temperature T = T c . Be- 
sides, when we select the value of the crystal field such as 
D/J = —1.38 and —1.4 we observe two successive phase 
transitions. In other words, if we cool the system starting 
from a finite temperature T > T c , the system undergoes 
a phase transition from paramagnetic to ferromagnetic 
phase at T — T c . If we keep on cooling process then 
the second order transition at a finite temperature is fol- 
lowed by a first order transition at a lower temperature 
T < T c . These results show the existence of reentrant 
phenomena. Furthermore, as the value of crystal field 
D/J decreases then the second order transition temper- 
ature decreases and the first order transition temperature 
increases. On the left panel in Fig. 2b, we see the be- 
havior of the transverse magnetization with temperature 
for some selected values of crystal field with h/J = 0. 
For D/J = 1.0,0,-0.5,-1.0 and -1.25 the transverse 
magnetization m x curves increase with the increase of 
the temperature and then show a cusp which increases 
in height as the value of the crystal field decreases at 
T = T c and decline as the temperature increases. Conse- 
quently, the transverse magnetization curves can be sepa- 
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FIG. 2: Temperature dependencies of (a) longitudinal mag- 
netization m z , (b) transverse magnetization m x , (c) longitu- 
dinal quadrupolar moment q z and (d) transverse quadrupolar 
moment q x for the spin-1 system on a honeycomb lattice with 
a fixed value of transverse field O/J = 0.5 and some selected 
values of crystal field D/J = f.0,0, -0.5, -1.0, -1.25, -1.38 
and —1.4. The longitudinal magnetic field value is selected 
as h/J = and 0.1 on the left and right panels, respectively. 
Dashed lines correspond to the first-order solutions. 



rated into two regions: the first is the nonmagnetic region 
in which m z = and the second is the magnetic region 
in which m z ^ 0. For D/J = —1.38 and —1.4 m x curves 
show a discontinuous behavior at the first order transi- 
tion point. In Figs. 2c and 2d, the variation of the longi- 
tudinal and transverse quadrupolar moments with tem- 
perature are shown, respectively. The same crystal field 
D/J values are used as in Figs. 2a and 2b. We see that 
the longitudinal quadrupolar moment q z decreases as the 
temperature increases and change abruptly at the second 
order transition temperature for D/J = 1.0, 0, —0.5, —1.0 
and —1.25. In case of D/J = —1.38 and —1.4, q z curves 
exhibit two minima which correspond to the first and 
second order transition temperatures. In contrast to q z 
curves in Fig. 2 c, the transverse quadrupolar moment q x 
curves increases as the temperature increases and change 
abruptly at the second order transition temperature for 
D/J = 1.0,0,-0.5,-1.0 and -1.25. For the values of 
crystal field D/J = —1.38 and —1.4 q x curves exhibit 
two maxima corresponding to the first and second order 
transitions in the system. When we apply a longitudinal 
magnetic field such as h/J = 0.1 on the system, both 
the first and second order phase transition effects on the 
order parameters m z , m Xl q z and q x are removed. This 
phenomenon can be clearly seen from the right hand side 
panels in Figs. 2a-2d. 

Next, the transverse field effect on the variation of the 
order parameters with temperature such as the longitudi- 
nal and transverse magnetizations m z and m x as well as 
the longitudinal and transverse quadrupolar moments q z 
and q x with D/J = —1.1 and h/J~0 can be seen on the 
left panels Fig. 3. Dashed lines represent the first order 
transition solutions. In Fig. 3a, we plot the longitudinal 
component of magnetization for selected values of trans- 
verse field fi/J = 0,0.5,0.75,0.9 and 1. As we can see 
on the left panel in Fig. 3a, as the transverse field fl/J 
increases, critical temperature value approaches to zero 
and saturation value of m z decreases. Hence, we can say 
that applying transverse field fi/J weakens the longitu- 
dinal component of magnetization. This is an expected 
result. However, two successive transitions (i.e. reentrant 
behavior) are observed on the longitudinal magnetization 
m z for Q/J = 0. It is clear that applied transverse field 
ft/ J destructs the first order transition. On the left panel 
in Fig. 3b, we represent the effect of the transverse field 
il/ J on the temperature dependence of m x . In contrast 
to the situation in m Zl transverse field f2/J strength- 
ens the transverse component of magnetization. We note 
that, at zero transverse field the transverse magnetization 
m x is zero for the whole range of temperature. On the 
left panels in Figs. 3c and 3d, we present the effect of the 
transverse field 0/ J on the variation of the quadrupolar 
moments q z and q x with temperature, respectively. As we 
can clearly see from these plots, longitudinal quadrupolar 
moment q z has two minima while transverse counterpart 
q x has two maxima corresponding to the first and second 
order transitions for tt/J = 0. Furthermore, applying 
any longitudinal magnetic field such as h/ J = 0.1 on the 
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0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 

kJ/J kJ/J 



(b) 




k ,™ k .™ 
(c) 




0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 



kJ/J kJ/J 

(d) 

FIG. 3: Temperature dependencies of (a) longitudinal mag- 
netization m z , (b) transverse magnetization m x , (c) longitu- 
dinal quadrupolar moment q z and (d) transverse quadrupo- 
lar moment q x for the spin-1 system on a honeycomb lattice 
with a fixed value of crystal field D/J = —1.1 and some se- 
lected values of transverse field Q/J = 0,0.5,0.75,0.9 and 
1.0. Q/J — 0.25 is also selected in (b). The longitudinal 
magnetic field value is selected as h/J = and 0.1 on the left 
and right panels, respectively. Dashed lines correspond to the 
first-order solutions. 



system, both the first and second order phase transitions 
are destructed. This behavior is illustrated on the right 
panels in Fig. 3. Hence, on the basis of these results 
(see the right panels in Figs. 2 and 3) we believe that 
the effects of the transverse field Q/J are very different 
from those of the longitudinal counterpart h/J because 
the origin of the transverse field is quantum mechanical 
and can produce quantum effects. 



Finally in Fig. 4, we present the numerical results 
for the temperature dependence of entropy per site and 
Helmholtz free energy of the system. As far as we know, 
there is not such a study dealing with the variation of the 
entropy of the system with the Hamiltonian parameters. 
In Figs. 4a and 4c, we show the effect of the crystal field 
D/J on the entropy and free energy for Q /J = and 
h/J = 0. For D/J = 0.5,0,-0.5,-0.75 and -1.0 en- 
tropy is not important at low temperatures and free en- 
ergy is equal to ground state energy of the system. How- 
ever, as the temperature increases, the system wants to 
maximize its entropy in order to minimize its free energy 
and hence entropy becomes important. Furthermore, the 
entropy of the system is continuous at the critical tem- 
perature which means that the type of the transition is 
second order. On the other hand, for D/J < —1.0 the 
entropy and free energy curves represented with dashed 
lines show a discontinuous behavior at low temperatures 
which indicates a first order transition and it originates 
from a discontinuous change in the internal energy of the 
system. Discontinuous behavior of the entropy can be 
clearly seen in the inset figure in Fig. 4a. As the absolute 
value of the crystal field increases, energy gap in free en- 
ergy curves (see Fig. 4c) disappear and ground state en- 
ergy equals to zero for tricritical crystal field value D t / J. 
Furthermore, we see that increasing the absolute value of 
the crystal field D/J makes the absolute value of free en- 
ergy decrease and causes the system to be in a consider- 
ably disordered state. These remarkable observations are 
not reported in the literature. The effect of the transverse 
field Q/ J on the temperature dependence of entropy and 
Helmholtz free energy can be seen in Figs. 4b and 4d 
for selected values of Q/J = 0,0.5,0.75,0.9 and 1.0 with 
D/J = —1.1 and h/J = 0. As the transverse field Q/J 
value increases, critical temperature decreases and the 
system reaches disordered phase and maximizes its en- 
tropy while minimizing its free energy at early stages of 
temperature scale. Besides, applied transverse field Q/J 
tends to destroy the energy gap in free energy and hence 
the first order transitions are removed. One should no- 
tice that although we observe the first order transition 
effects for Q/J = in free energy (see Fig. 4c), this 
effect is not evident for the entropy (Fig. 4b). This is 
due to the fact that, when we select D/J = —1.1 and 
h/J = with Q/J — 0, the system is in a disordered 
state and the internal energy is zero at low temperatures 
(i.e. ground state energy) and its value is not changed 
until the system undergoes a first order transition. 
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IV. CONCLUSION 




(a) 




(b) 




(c) 




In this work, we have studied the phase diagrams of the 
spin-1 transverse Ising model with a longitudinal crystal 
field in the presence of a longitudinal magnetic field on 
a honeycomb lattice within the framework of the IEFT 
approximation that takes into account the correlations 
between different spins in the cluster of considered lattice. 
We have given the proper phase diagrams, especially the 
first-order transition lines that include reentrant phase 
transition regions. Both the order parameters (m z , m x , 
Qzi Qx) and Hclmholtz free energy F and the entropy 
S curves show discontinuous and unstable points of the 
system and support the predictions in Figs, la and lb. 

A number of interesting phenomena such as reentrant 
phenomena have been found in the physical quantities 
originating from the crystal field as well as the trans- 
verse and longitudinal components of the magnetic field. 
We have found that one can observe reentrant behav- 
ior in the system for the values of 51/ J < 0.861 and 
— 1.4533 < D/J < —1.0201 and the tricritical points ex- 
ist for 1.4333 < -D t /J < 1.2776 and fi/J < 0.861. The 
results show that the reentrant phenomenon originates 
from the competition between the crystal field D/J and 
transverse field ft/ J. Besides, applying a transverse field 
Q/ J on the system has the tendency to destruct the first 
order transitions, while the longitudinal counterpart h/J 
destructs both the first and second order phase transi- 
tions. Hence, we believe that the effects of the transverse 
field fi/J are very different from those of the longitudi- 
nal counterpart h/J since the transverse field can pro- 
duce quantum effects. These interesting results are not 
reported in the literature. 

We hope that the results obtained in this work may be 
beneficial from both theoretical and experimental point 
of view. 
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(d) 



FIG. 4: Entropy of the system as a function of temperature 
for h/J = and (a) fi/J = with some selected values of 
crystal field D/J = 0.5, 0, -0.5, -0.75, -1.0, -1.4 and -1.44, 
(b) fi/J = 0, 0.5, 0.75, 0.9 and 1.0 with a fixed value of crystal 
field D/J = —1.1. (c) Temperature dependence of Helmholtz 
free energy of the system for h/ J = 0. The value of transverse 
field is fixed as fi/J = 0.0 with some selected values of crystal 
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Appendix A 

The basis relations corresponding to equations dHJ) , (O , 
1(12)1 CEU) (H5]) and (HHJ). The coefficients k, p i: r u v t 
(i = 0,1,..., 5); kj, Cj, rij, fXj (j = 1,...,4); and a k , b k 
(k = 1, 2, 3) can be derived from a mathematical identity 
exp(aV)F(x) = F(x + a). 

(5 z ) = i +3/ci<Si)+3(Ji-io)<S?)+3i 2 (Si5 2 ) 

+6(fe 2 - fci)(Si^ 2 2 ) + 3{l - 2h + l 3 )(SfSl) 
+h(S 1 S 2 S 3 ) + 3(Z 4 - l2)(S!S 2 S 2 3 ) 
+3(k 1 ~2k 2 + k 4 )(S 1 SiS 2 ) 
+(-l + 3/i - 3Z 3 + l 5 )(S 2 S 2 Sl) (Al) 

h = F(0) 

h = cosh(JV)F( X ) \x=0 

h = sinh 2 (JV)F(x)| a;= o 

h = cosh 2 {JV)F(x)\ x=0 

U = cosh(JV)smh 2 (JV)F(a;)| x= o 

l 5 = cosh 3 (JV)F(x)| x=0 

fci = sinh(JV)F(x)| a;= o 

fc 2 = cosh(JV)sinh(JV)F(x)| a:= o 

fc 3 = sinh 3 (JV)F(a;)| x= o 

fc 4 = cosh 2 (JV)sinh(JV) J F(x)| :c= o 



(So) = Po + 3c 1 (S 1 )+3(p 1 -p )(S 2 )+3p 2 (S 1 S 2 ) 
+6(c 2 - ci)(SiS|) + 3(p - 2 P i + p 3 )(S 2 Sf) 
+c 3 (SiS 2 S 3 ) + 3(p 4 -p 2 )(SiS 2 S 2 ) 
+3(c 1 -2c 2 + c 4 )(Si,S' 2 2 S 2 ) 
+(-Po + 3pi - 3 P3 + P5 )(S 1 2 S 2 ^3 2 ) (A2) 

Po = H(0) 

Pi = cosh(JV)i2"(a;)| x =o 

p 2 = smh 2 (J\7)H(x)\ x=0 

p 3 = cosh 2 (JV)i? (a;) | x=0 

Pi = cosh(JV)sinh 2 ( JV)i? (a;) | x=0 

p 5 = cosh 3 (JW)H(x) \ x=0 

d = sinh(JV)ff(a;)| x= o 

c 2 = cosh( JV)sinh( JV)H(x)\ x= o 

c 3 = smh a (JV)H(x)\ x =o 

c 4 = cosh 2 (JV)sinh(JV)i/( X )\x=0 

(SO = ai (1 - ((S 2 ) 2 )) + a 2 (S%) + a 3 ((S*) 2 ) (A3) 



01 = F( 7 ) 

a 2 = sinh(JV)F(x + 7)U=o 
a 3 = cosh(JV)F(a; + 7)| x=0 

((S 2 ) 2 ) = r + 3n 1 (5 1 )+3(r 1 -r )(5 1 2 )+3r 2 (5 1 5 2 ) 

+6(n 2 - ri^Sf) + 3(r - 2r x + r 3 )(S 2 S 2 2 ) 
+n 3 (S'iS' 2 S'3) + 3(r 4 - r 2 ){S 1 S 2 S 2 ) 
+3{n 1 -2n 2 + n i )(S 1 S 2 2 Si) 
+ (-r + 3ri - 3r 3 + r 5 )(S 2 5 2 ^3 2 ) (A4) 

r = G(0) 

r x = cosh(JV)G(a;)| x =o 

r 2 = sinh 2 (JV)G(a;)| K= o 

r 3 = cosh 2 (JV)G(a;)U = o 

r 4 = cosh(JV)sinh 2 (JV)G(a;)| x =o 

r s = cosh 3 (JV)G(jc) | x =o 

ni = sinh( JV)G(x)| x= o 

n 2 = cosh(JV)sinh(JV)G(a;)| x= o 

n 3 = sinh 3 (JV)G(a;)| x= o 

n 4 = cosh 2 ( JV)sinh( JV)G(a:)| a;= o 

((^) 2 ) = i;o + 3/ Jl (5 1 )+3(7 Jl - Wo )(^ 1 2 )+3« 2 (5 1 52) 

+6(^ 2 - Mi)(S!5 2 ) + 3{v - 2«i + v 3 )(S 2 S 2 ) 
+p, 3 (SiS 2 S 3 ) + 3(u 4 - v 2 )(S 1 S 2 S 3 1 } 
+3{p 1 -2p 2 + p 4 )(S 1 S 2 2 S 2 ) 
+ (-v a + 3tji - 3t> 3 + v 5 )(S 2 S 2 S 2 ) (A5) 





- K(0) 


Vl 


= cosh(JV)if(a;)| 2;= o 


V2 


= sinh 2 (JV)if(a;)U = o 


V3 


= cosh 2 (JV)^(a;)U = o 


Vi 


= cosh(JV)sinh 2 (JV)X(x)| : 


V5 


= cosh 3 (JV)if(a;)U = o 


Ml 


= smh{J\7)K(x)\ x = 


M2 


= cosh(JV)sinh(JV)iT(a;)| x: 


M3 


= sinh 3 (JV)X(a;)U = o 


/'4 


= cosh 2 (JV)sinh(JV)/sT(a;)| : 



(S 2 ) = h(l- ((^) 2 )) + b 2 (SZ) + 6 3 ((S 2 ) 2 ) (A6) 
h = G( 7 ) 

b 2 = sinh(JV)G(a; + 7)| x= o 
63 = cosh(JV)G(ir + 7)| x= o 
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The complete set of twenty three linear equations of the 
spin-1 honeycomb lattice (q = 3): 



(S%) = i + 3fci<Si) 



(Si So} 



(S1S2S0) 



(Si) 
(SiS 2 ) 
(Sls 2 s 3 ) 

(SI) 
(SiS 2 2 ) 

(sis!) 

(S Sf) 
(SoS\S 2 ) 
(SqSiS 2 ) 



(S\S 2 S 3 ) 
(SiS 2 S 3 ) 

(sfslst) 

((Stf) 



+6(fc 2 -£ 4 )<S 1 S 2 2 } 



3(h-l )(S 2 ) + 3l 2 (S 1 S 2 ) 
+ 3(l -2h + l 3 )(S 2 1 Sl) 
+fc 3 (SiS 2 S 3 ) + 3(Z 4 - ^(S^Sl) 
+3(k 1 -2k 2 + k i )(S 1 SiSl) 

+(-i + 3h-3i 3 + k)(s 2 s 2 s 2 ) 

{3l 1 -2l )(S 1 ) + 3k 1 (Sf) 
+3(l -2l 1 + l 2 + l 3 )(S 1 S 2 ) 
+6(fc 2 - k^SfSl) + fc 3 <Si5 2 S|) 
+Ho + 3ii - 3Z 2 - 3/3 + 3Z 4 + Ib)(SiS%S£) 
+3(k 1 -2k 2 + k 4 )(S 2 1 SlSl) 
(Z - 3/i + 3Z 2 + 3i 3 )<SiS 2 ) + (6fc 2 - 3fci)(Si5 2 ) 
+(-/ + 3/i - 3Z 2 - 3l 3 + 3l 4 + k)(SiS 2 Sl) 
+(3fci - 6fc 2 + fc 3 + 3k 4 )(S 1 S 2 1 S 3 2 ) 
ai (l-((5 z ) 2 ))+a 2 (^)+a 3 ((5 z ) 2 ) 
a 1 (5 1 )+a 2 (5o5 1 ) + (a 3 -a 1 )(5 1 5 2 > 
ai <SiS 2 ) + a 2 (S Q S 1 S 2 ) + (a 3 - a 1 )(5 1 5 2 5 2 ) 
fe 1 (l-((5 z ) 2 » + fe 2 (5 ( 5)+6 3 ((5 z ) 2 ) 
&i(5i)+6 2 (5 5 1 ) + (fe 3 -6 1 )(5 1 5 2 > 
&i(5 2 ) + 6 2 <S'oSi> + (63 - h)(S 2 S 2 ) 



<SiS|Sg> = (r -3ri+3r2 + 3r 3 )<SiS|) 
+ (-3n 1 +6n 2 )(5 1 5 2 ) 

(-r + 3n - 3r 2 - 3r 3 + 3r 4 + r 5 )(SiS 2 2 S 3 2 ) 
+(3m - 6n 2 + n 3 + 3n 4 )(S 1 S 2 S 3 l ) 
(S 2 S 2 S 2 ) = (ro-3r 1 + 3r 2 + 3r 3 )(S 2 S 2 ) 
+(-3ni+6n 2 )<SiS|> 
+ (3ni — 6n 2 + n 3 + 3ri4)(S'i5|S'|) 
+(-r + 3n - 3r 2 - 3r 3 + 3r 4 + r 5 )(S 2 S 2 S 2 ) 
(SS) = Po+3c 1 (S 1 )+3(p 1 -p Q )(S 2 1 ) + 3 P2 (S V S 2 ) 
+6(c 2 - ci)(5 1 5 2 > + 3(p Q - 2 Pl + p 3 )(S 2 S 2 ) 
+c 3 (S 1 S 2 S 3 ) + 3(p 4 -P2) (SiS 2 S 2 ) 
+3(c 1 -2c 2 +c 4 )(S 1 S 2 S* 2 ) 
+(-Po + 3pi - 3p 3 +p 5 )(5 1 2 5|5|) 
((S£) 2 ) - Wo + 3 A11 (^ 1 >+3( Ul -«o)(5 1 2 )+3 U2 (5 1 5 2 ) 

+6(^ 2 - ^i)(Si5 2 > + 3(w - 2ui + v 3 )(S 2 S 2 ) 
+H 3 (S 1 S 2 S 3 ) + 3(v 4 - v 2 )(SiS 2 S 3 ) 



h)(S 1 S 2 S 2 ) 
■6i)<Si5|Sg) 



& 3 (5 )+6 2 (5 2 ) 

& 3 (SoSi)+& 2 (SiS 2 ) 
& 3 (5o5 2 )+6 2 (5 2 5 2 ) 
&i(5i^ 2 )+6 2 (5 5 1 5 2 ) + (fe 3 - 

&i(5 2 5|) + b 2 (S S 2 S 2 ) + (63 - &i)(S 2 S 2 S 2 ) 
r + 3m<Si) + 3(n - r )<S?) + 3r 2 <SiS 2 ) 
+6(n 2 - m)<SiSf ) + 3(r - 2n + r 3 )(S 1 2 S 2 2 ) 
+n 3 ((^5 2 5 3 ) + 3(r 4 - r 2 )(SiS 2 Sf ) 
+3(n 1 -2n 2 + n 4 )(5i5 2 2 5* 2 ) 
+(-ro + 3r 1 -3r 3 + r 5 )(S 1 2 S 2 S 2 ) 
(SiSg) = (3r 1 -2r )(5 1 )+3n 1 (5 1 2 ) 

+(3r 2 + 3r -6r 1 +3r 3 )(5i5 2 ) 
+6(n 2 - ni )(S 2 S 2 ) + n 3 <SiS 2 S3> 
+(-r + 3n - 3r 2 - 3r 3 + 3r 4 + r 5 )(SiS 2 S 3 2 ) 
+3(n 1 -2n 2 + n 4 )(5 1 2 S*|5|) 
(S?Sg) - (3r 1 -2r )(5 1 2 )+3n 1 (5 1 ) 

+(3r 2 + 3r -6r 1 +3r 3 )(S 1 2 S 2 ) 
+6(n 2 -n 1 )(S 1 S 2 ) 
+(3ni — 6n 2 + n 3 + 3n 4 )(5i5|S'|) 
+(-r + 3n - 3r 2 - 3r 3 + 3r 4 + r B )<S?Sf Sf) 
(S V S 2 S 2 ) = (ro-3ri+3r 2 + 3r 3 )<SiS2> 
+(-3ni + 6n 2 )(5i5 2 ) 

(— r + 3ri — 3r 2 - 3r 3 + 3r 4 + r 5 )(SiS 2 S 3 ) 
+(3ni — 6n 2 + n 3 + 3n 4 )(Si5f S* 2 ) 



+3(/ii - 2^ 2 
+(-i>o + 3ui 



H 4 )(SiS 2 S 2 



3v 3 



v 5 )(S 2 S 2 u 3 



si) 



(A7) 
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